#============================================================================
#Modeling the transmission of Chagas Disease between triatomines and cuy.
#Built from the skeleton up, starting 08/13/10
# pushed to Google Code on Feb 14th, valentine's day, 2011
# such a nice spring day in the middle of the winter
#============================================================================
  
  #==========================================================
  #Import CSV with run parameters
  #==========================================================
  inputdir<-"/Users/mzlevy/Documents/Superspreading/"
  outfolder<-"feb14_alt10"
  outdir<-paste(inputdir,outfolder, sep="")
  setwd(inputdir)  #set working directory to retrieve data
  pathdata<-"Run_Tracking_feb14_mzl_alt10.csv"
  runTracking<-read.csv(pathdata)     #Import csv with run parameters
  print(dir.create(paste(inputdir,outfolder, sep=""), recursive=TRUE))
  
  
  #==========================================================
  #Assign columns from csv to vectors
  #==========================================================
  VrunNumber<<-runTracking$RunNumber	       #Import vector of run numbers
  VNbugs<<-runTracking$Nbugs	               #Import vector of number of bugs
  VNcuy<<-runTracking$Ncuy 	                 #Import vector of number of guinea pigs
  VPercentBugI<<-runTracking$InitiallyInfectedBugs 	 #Import vector of the percent of bugs initially infected
  VPercentCuyI<<-runTracking$InitiallyInfectedCuy	   #Import vector of the percent of cuy initially infected
  #VBNlength<<-runTracking$BottleneckLength	 #Import vector of bottleneck lengths
  #VBNstart<<-runTracking$BottleneckStart	   #Import vector of bottleneck start days
  VTranCtoB<<-runTracking$TransmissionCtoB	 #Import vector of transmission from cuy to bug rates
  VTranBtoC<<-runTracking$TransmissionBtoC	 #Import vector of trasmission from buy to cuy rates
  VcLife<<-runTracking$LifespanCuy	         #Import vector of the number of days cuy live
  VbLife<<-runTracking$LifespanBugs          #Import vector of the number of days bugs live

#==========================================================
#Set initial parameters to the parameters from CSV in row p
#==========================================================

for (p in 1:length(VrunNumber))   
{                                   

#==============================
#Parameters sets all the initial parameters as globals
#==============================

parameters<-function()
{
#==========================
#Initial parameters
#==========================
  Ncuy<<-VNcuy[p]                                     #Number of cuy
  Nbug<<-VNbugs[p]                                    #Number of bugs
  NcuyI<<-(VNcuy[p]*VPercentCuyI[p])                  #Number of infected cuy
  NbugI<<-(VNbugs[p]*VPercentBugI[p])                 #Number of infected bugs

#===========================
#Parameters for when bugs feed
#===========================
  f<<-40              #Number of times bugs feed, columns of the FeedDay matrix,
  l<<-30              #Number of days between bloodmeals for each bug
  scalel<<-.005       #Bigger means more variance
  shapel<<-l/scalel

#=============================
#Parameters for the frequency with which vectors and hosts die
#=============================
  cLife<<-VcLife[p]                       #frequency with which host are reset
  bLife<<-VbLife[p]                      #frequency with which vectors are reset
  
#=============================
#Monod function and parameters for BugtoCuy infectiousness
#=============================
  da<-seq(1, l*f)
  r<<-VTranBtoC[p]                          #r is the maximum infectiousness
  a<<- 9                                    #a is the # days to get up to half as infectious as the max
  InfectProbBugs<<-((r*da)/(a+da))     #Vector of bug infectiousness 
  
#==============================
#Creates a vector of solutions for cuy infectiousness.  The solution run
#between day 1 and the last day of the simulation, day l*f.  InfectProb1 will be
#used for <=20 days post infection and InfectProb2 will be use for >20 days
#post infection.  This is a simple fit to data.
#==============================
  InfectProb1<<-(VTranCtoB[p]*exp(0.33*(da[1:20]+10))/exp(0.33*(20+10)))  #Infectiousness for days 1-20 post-infection
  InfectProb2<<-(VTranCtoB[p]*exp(-0.0603*(da[21:length(da)]-180))/exp(0.33*(20+10))) #Infectiousness for days 21 on post-infection
  InfectProbCuy<<-c(InfectProb1, InfectProb2)   #Vector for cuy infectiousness 
  change<-which(InfectProbCuy>1)                
  InfectProbCuy[change]<<-1                     #Changes values over 1 to 1
      
#==============================
#pairfunction parameters for bottlenecks
#==============================
  #startBottle<<-VBNstart[p]                 #Day to start the bottleneck
  #lengthBottle<<-VBNlength[p]               #Number of days in in the bottleneck
  #endBottle<<-startBottle+lengthBottle      #Day to end the bottleneck 
  #bottleNeck<<-2                            #Number of cuy to stay active during a bottleneck
  EndTime<<-l*f
  #notInBottle<-rep(1, Ncuy)
  #notInBottle<<-notInBottle(-hostBottle)
  
#=============================
#Create the matrices that will hold information regarding 
#which bug infected which cuy (and vice versa) and on what day
#=============================
  bugToCuy<<-matrix(0, nrow=0, ncol=3)     #col 1= bug number, col 2= cuy number, col 3= day
  cuyToBug<<-matrix(0, nrow=0, ncol=3)     #col 1= cuy number, col 2= bug number, col 3= day
  

} #end parameters()

#==========================================================================
#initialize is a function that is called at the start of every run to set the 
#run specific parameters 
#==========================================================================

initialize<-function()
{

#============================
#Sets the initial parameters for the host infection status, cI, host infection 
#day, cID, which host are active/ alive, activeHost and the total number of 
#active hosts that are infected.
#============================
    cI<<-rep(0, Ncuy)     #create a 0/1 vector to record host infection
    cID<<-rep(0, Ncuy)    #create a vector to record host infection day
    activeHost<<-seq(1, Ncuy) #create a vector of avaliable/ living hosts
    cuyTime<<-rep(0, EndTime)       #A matrix to hold the daily sum of cuy infected
    phi<-sample(1:Ncuy, NcuyI, replace=FALSE)       #Which cuy to be infected
    cID[phi]<<-1                                    #Assign the infection day to day 1
    cI[phi]<<-1


#============================
#Sets the initial parameters for the vector infection status, bI, the vector
#infection day, bID, which vectors are active/ alive  and the total number of 
#active vectors that are infected on a given day
#============================
    bI<<-rep(0, Nbug)
    bID<<-rep(0, Nbug)
    activeVector<<-seq(1, Nbug)
    bugTime<<-rep(0, EndTime)       #A matrix to hold the daily sum of bugs infected
    phi<-sample(1:Nbug, NbugI , replace=FALSE)      #Which bug to be infected
    bID[phi]<<-1                                    #Assign the infection day to day 1
    bI[phi]<<-1

#================================
#Initialize calls the bugDeath and cuyDeath functions which create the 
#matrices of assigned death days for each host and vector
#================================
    bDeath<<-bugDeath()
    cDeath<<-cuyDeath()

} #end initialize()

#==================================
#Creates a matrix that pre-determines the death day of each vector.  Death 
#frequency is determined by bLife. This function returns a matrix bDeath with a 
#1 assigned to the bug's death day.  
#==================================

bugDeath<-function()
    {
    deathFreq<-1:ceiling((l*f)/bLife)*bLife  #Determines the number of times each bug could die considering 
                                             #the length of a bug's life, bLife and the number of days in the simulation, l*f. 
                                             #It then produces a vector of the last possible death day for each death cycle 
    bDeath<-matrix(0, nrow=Nbug, ncol=max(deathFreq)) #Matrix with columns = days and vectors = rows
    for (k in 1:Nbug)
      {
      eL<-deathFreq-sample(0:(bLife-1),1) #creates a vector of days vector k will die
      bDeath[k, eL]<-1                  #assigns vector of days to bDeath matrix
      } #end outter for loop
      return(bDeath)
    } #end bugDeath

#==================================
#cuyDeath is a function that creates a matrix that pre-determines the death day 
#of each host.  The start day is random but then they are reset every cLife 
#number of days therefore cLife controls the frequency of death.  
#==================================

cuyDeath<-function()
  {
  deathFreq<-1:ceiling((l*f)/cLife)*cLife
  cDeath<-matrix(0, nrow=Ncuy, ncol=max(deathFreq))
  for (k in 1:Ncuy)
    {
      eL<-deathFreq-sample(0:(cLife-1),1) #creates a vector of when hosts will die
      cDeath[k, eL]<-1   #Assigns vector of death days to cDeath matrix
    } #end outter for loop
    return(cDeath)
  } #end cuyDeath()
  
 #=================================
 #resetBug is a function that looks at column day of the cDeath matrix to see which 
 #bugs are due to be reset on that day and resets bI and bID to 0.
 #=================================
      resetBugs<-function(day)
        {
        if (sum(bDeath[,day])==0) 
          {
          return          #Return if no bugs die that day
          } #end if
        else
          {
            reset<-which(bDeath[,day]==1) #Rows in bDeath of vectors to die that day
            newVector<-seq(length(bI)+1, length(reset)+length(bI)) #Creates a sequence of the new vectors 
            activeVector[reset]<<-newVector    #Replaces vectors that died with new vectors
            bID<<-c(bID, rep(0, length(reset)))#Adds new vectors to the end of bug infection day, bID
            bI<<-c(bI, rep(0, length(reset))) #Adds new vectors to the end of bug infection status, bI
          } #end else
        } #end resetBugs()
 #=================================
 #resetCuy is a function that looks at column day of the cDeath matrix to see which
 #cuy are due to be reset on that day and resets cI and cID to 0.
 #=================================
      resetCuy<-function(day)
        {
        if (sum(cDeath[,day])==0) 
          {
          return
          } #end if
        else
          {
            reset<-which(cDeath[,day]==1) #Rows in cDeath of hosts to die that day
            newHost<-seq(length(cI)+1, length(reset)+length(cI))  #Creates a sequence of the new hosts
            activeHost[reset]<<-newHost #Replaces hosts that died with new hosts
            cID<<-c(cID, rep(0, length(reset)))#Adds new hosts to the end of cuy infection day, cID
            cI<<-c(cI, rep(0, length(reset))) #Adds new hosts to the end of cuy infection status, cI
          } #end else
        } #end resetBuy()
 #=============================
 #Function whichFeed determines which bugs are feeding on day i
 #whichFeed returns FT, a vector of bug numbers or NA if no bugs are feeding
 #=============================
      whichFeed<-function(Feed, day)
        {
        total<-sum(Feed[,day])
        if (total==0) return(NA)
        rowsFeeding<-which(Feed[,day]==1)#determines the rows in column day that have a 1
        FT<-activeVector[rowsFeeding] #determines the vectors that coorispond to the rows
        return(FT)  #returns a vector of vectors feeding that day
        } #end whichFeed
 #=============================
 #Function pairfunction matches bugs that are feeding on a given day day 
 #(vectors in the vector FT) with potential host.  pairfunction returns a matrix
 #with column 1= vectors and  #column 2= the host that vector will feed on.
 #Need to input pairtype as "bottleNeck" or "random" 
 #=============================
      pairfunction<-function(FT, day, pairtype)
        {
          BugHost<-matrix(NA, nrow=length(FT), ncol=2) #Create matrix of vectors (col 1) and hosts (col 2)
          BugHost[,1]<-FT                              #Assign vectors feeding on day to column 1
   #==================================
   #No bottleneck or special pairing option
   #all pairing is random and even
   #==================================
     if (pairtype=="random")
     {
        BugHost[,2]<-sample(activeHost, length(FT), replace=TRUE)  #Assigns cuy to the second column
        return(BugHost)
     } #end if for pairtype "random"
   #==================================
   #Bottleneck option
   #===================================
     if (pairtype=="bottle")
     {   
        if (day == endBottle)
          {
            resumeHost()#On the last day of the bottleneck, cI and cID are reset
          } #end if
        if (day == startBottle)
          {
            hostBottle<<-sample(activeHost, bottleNeck, replace=FALSE)    #determines which cuy will remain "active"  
            BugHost[,2]<-sample(hostBottle,length(FT), replace=TRUE)
            return(BugHost)
          }
        if ((day <= startBottle) | (day >= endBottle))         #Gives time constraints to bottleneck
        {
          BugHost[,2]<-sample(activeHost,length(FT), replace=TRUE)  #Assigns host to the second column randomly if not in the bottleneck
          return(BugHost)
        } #end if
        else
          {
            BugHost[,2]<-sample(hostBottle,length(FT), replace=TRUE)
            return(BugHost)
          } #end else
     } #end if for pairtype "bottle"
  } #pairfuntion()	
  #================================================
  #resumeHost is a function that is called at the end of a bottleneck.
  #It puts the cuy infection day (cID) back to 0 for all cuy that were removed
  #and makes the cuy infection status (cI) 0.
  #================================================
  
  resumeHost<-function()
        {
        cleancI<-rep(0, length(cI))                       #place holder vector for cI
        cleancID<-rep(NA, length(cID))                     #place holder vector for cID
        for (w in 1:length(hostBottle))                 #Loops through the active cuy
          {
            cleancI[hostBottle[w]]<-cI[hostBottle[w]]       #Keeps cI status for active cuy
            if (is.na(cID[hostBottle[w]])==FALSE)       #Keeps cID for active cuy
              {
                cleancID[hostBottle[w]]<-cID[hostBottle[w]]
              } #end if
          } #end for loop
            cI<<-cleancI                            #Resets cI
            cID<<-cleancID                          #Resets cID
        } #end zeroCuy	
	
 #=============================
 #transmissionSort is a loop function that takes in BugHost and day.  It 
 #determines if the vector/host pair should be passed to the InfectBtoC, 
 #InfectCtoB or neither
 #=============================

    transmissionSort<-function(BugHost, day)
      {
      DAY<<-day
      for(alpha in 1:length(BugHost[,1]))
        {
        
        vector<- BugHost[alpha, 1]
        host<- BugHost[alpha, 2]
        
        if ((bI[vector] == 1) & (cI[host] == 0) & (DAY>bID[vector]))#Vector that has been infected for at least one day an uninfected host
          {
            InfectBtoC(vector, host, day)      ##Pass the vector, host and day on to InfectBtoC()           
          } #end if
        if((bI[vector]==0) & (cI[host]==1) & (DAY>cID[host]))#Host that has been infected for at least one day and an uninfected vector
          {                                              
            InfectCtoB(vector, host, day)      #Pass the vector, host and day on to InfectCtoB()
          } #end if
        } #end for
      } #end preInfect()

#============================================================================
#InfectBtoC is a function that determines if infection is passed from an 
#infected vector to a non-infected host.  It takes in the day, vector number and 
#host number.  The function determines if the infection was passed between 
#vector and host.  If so, cI, cID bugToCuy are altered. 
#============================================================================

    InfectBtoC<-function(vector, host, day)
      {
      
      daysSinceInfection<-(day-bID[vector])  #Days since vector was initially infected
      InfectProb<-InfectProbBugs[daysSinceInfection]  #InfectProbBugs is a vector of solutions to the bug infectiousness equations

      randomNum<-runif(1, min=0, max=1)         #Generates a random number 
      if (randomNum <= InfectProb)              #If the random number is < the bug's infectiousness, the cuy is infected
        {
        cID[host]<<- day                  #Cuy infection day, cID is set
        cI[host]<<- 1                     #Cuys is set as infected
        bugToCuy<<-rbind(bugToCuy, c(vector, host, day)) #Transmission is added to the tree
        } #end if
      } #end InfectBtoC

#==============================
#InfectCtoB is a function that determines if infection is passed from an 
#infected host to a non-infected vector.  It takes in the day, vector number and 
#host number.  The function determines if the infection was passed between 
#vector and host.  If so, bI, bID cuyToBug are altered. 
#==============================

    InfectCtoB<-function(vector, host, day)
      {

          daysSinceInfection<-(day-cID[host])  #Days since host was initially infected
          randomNum<-runif(1, min=0, max=1)    #Generates a random number
          
          if(randomNum <= InfectProbCuy[daysSinceInfection])#InfectProbCuy is a vector of solutions to the host infectiousness equation
            {
            bID[vector]<<-day  #If the vector becomes infected, that vector's infection day is altered
            bI[vector]<<-1     #The vector's infection status is chaged from 0 to 1
            cuyToBug<<-rbind(cuyToBug, c(host, vector, day)) #The transmission event is added to the tree
            } #end if
      } #end InfectCtoB


#==========================
# setFeed is a function taht creates a matrix "Feed".  The days each vector will
#feed are pre-determined and indicated in the matrix Feed with a 1.
#feedtype can be "random", "sync" or "gamma"
#==============================

setFeed<-function(feedtype)
    {
    Feed<-matrix(0,nrow=Nbug, ncol=l*f)  #set Feed matrix
    everyl<-1:f*l                        #Make even every l days matrix
    #===================================================
    #create Feed for feedtype = "random"
    #===================================================
    if (feedtype=="random")
        {

            for (k in 1:Nbug)            
            {
                eL<-everyl-sample(0:(l-1),1)
                for (j in 1:f)
                {
                    Feed[k,eL[j]]<-1
                } #end interior for loop
            } #end exterior for loop
       return(Feed)
        } #end if feedtype = random
    #===================================================
    #create Feed for feedtype = "sync"
    #===================================================
    if (feedtype=="sync")
        {

            for (k in 1:Nbug)            
            {
                everyl<-(1:f*l)          #Creates a vector of what day ALL the vectos will feed
                Feed[,everyl]<-1         #Puts a 1 in the column of Feed when the vecotrs will be feeding
            } #end for loop
        return(Feed)
        } #end feedtype = sync
    #===================================================
    #create Feed for feedtype = "gamma"
    #===================================================
    if (feedtype=="gamma")
      {
        Feed<-StartFeed()  #StartFeed() is only used to create the feed function for a gamma feed pattern
        return(Feed)
      } #end if feedtype = gamma
      
    print("bad feed name") #Catch incase feedtype is invalid
    } #end setFeed() 


#============================
#StartFeed creates a feed day matrix Feed using a gamma distribution.
#StartFeed is called in the setFeed("gamma") function
#============================

StartFeed<-function()
{

    FDInt<-round(matrix(c(rgamma(1:1000, shape=shapel, scale=scalel)), nrow=Nbug, ncol=f))
                                                  #Create a matrix with feed day intervals (rounded to integers)
    LastFD<-matrix(nrow=Nbug, ncol=f)             #Create LastFeedDay matrix
    RowSum<-rowSums(FDInt)
    LastFD[,f]<-RowSum                            #Sum of the intervals is the last time each bug will feed
    omega<-which.max(RowSum)
    Ncolumns<-RowSum[omega]
    p <- (f-1)                                    #Start with column on to the left of the final column
    while ( p > 0 )
    {                                             
        LastFD[,p]<- RowSum-FDInt[,(p+1)]         #Assigns feed days for each vector
        RowSum<-LastFD[,p]

        p <-(p-1)
    } #end while

    Feed<-matrix(0, nrow=Nbug, ncol=Ncolumns)     #Matrix of 1/0 days each bug fed

    for (m in 1:Nbug)                             #Fills in the binary matrix Feed for days each bug feeds
    {                                             #the loop goes through rows first then columns
        for(n in 1:f)
        {
            LastFD[m,n]->x
            Feed[m,x]<-1                   
        } #end interior for loop
    } #end exterior for loop
return(Feed)
} #end StartFeed()

#============================
#Main loop
#============================
runall<-function(Feed)
{
    #initialize()
    for (day in 1:EndTime)
        {
        FT<- whichFeed(Feed, day)               #Returns bugs feeding on day i as a vector of bug numbers
        if(is.na(FT[1]) == TRUE)
          {
            return             #Goes to next day if no bugs are feeding on day i
          } #end if
            else
            {
                BugHost<<-pairfunction(FT, day, "random")       #Pairs bugs and hosts in a matrix, column 1= bugs, column 2= host
                                                             #Use "bottle" or "random"
                transmissionSort(BugHost, day)  #Determines if host-vector pair 
                                                #should be sent to InfectBtoC, 
                                                #InfectCtoB or neither. 
             } #end else                        
         cuyTime[day]<-sum(cI[activeHost])   #Records the number of hosts infected at the end of each day
         bugTime[day]<-sum(bI[activeVector])   #Records the number of vectors infected at the end of each day
         #cuyDay[1:length(cI),]<<-cI                  #A matrix of each cuy's infection status each day
         #bugDay[1:length(bI),]<<-bI                  #A matrix of each bug's infection status each day
         #for (y in 1:Ncuy)
         #   {
         #     if (day>=startBottle && day<=(startBottle + lengthBottle) && is.na(notActive[y])==TRUE)
         #       {                             #=============================
         #         cuyDay[day, y]<<-999        #Comment out when not running bottleneck
         #       } #end if                     #=============================
         #   } #end for              
                  
         resetBugs(day)
         resetCuy(day)
        } #end for
   
    return(cbind(bugTime,cuyTime))
} #end runall()

#================================
#Loop used to call the whole program.
#It first calls the function to set all of the parameters as globals.
#================================
nsims<-100                                 #number of simulations to run
parameters()
out<-matrix(NA, nrow=EndTime,ncol=2)                               #Sets the global parameters
outVector<-matrix(NA, nrow=nsims,ncol=EndTime)
outHost<-matrix(NA, nrow=nsims,ncol=(l*f))
maxcuy<-matrix(NA, nrow=nsims, ncol=2)  #Creates matrix that will hold the max number of hosts infected during a simulation
maxbugs<-matrix(NA, nrow=nsims, ncol=2) #Creates matrix that will hold the max number of vectors infected during a simulation
cIDmatrix<-matrix(nrow=nsims, ncol=Ncuy)#Creates matrix that will hold the infection day of each host during a simulation
bIDmatrix<-matrix(nrow=nsims, ncol=Nbug)#Creates matrix that will hold the infection day of each vector during a simulation
#cuyDay<<-matrix(NA, nrow=EndTime, ncol=sum(colSums(cDeath)))    #Infection status of each cuy each day
#bugDay<<-matrix(NA, nrow=EndTime, ncol=sum(colSums(bDeath)))    #Infection status of each bug each day
 
for (simulationNumber in 1:nsims)
    {
    initialize()                    #Calls and sets parameters from csv
    Feed<-setFeed("random")         #Sets feed type.  "random", "gamma", "sync"
    out<-runall(Feed)           #Saves the results of each simulation in out
    outVector[simulationNumber,]<-out[,1]
    outHost[simulationNumber,]<-out[,2]
    cIDmatrix[simulationNumber,]<-cID[activeHost] #Saves final cID for that simulation
    bIDmatrix[simulationNumber,]<-bID[activeVector] #Saves final vID for that simulation
    maxcuy[simulationNumber]<-max(cuyTime)   #Saves the max number of hosts infected over the simulation
    maxbugs[simulationNumber]<-max(bugTime)  #Saves the max number of vectors infected over the simulation
    print(simulationNumber)
    } #end for loop
    
#========================
#Plotting Commands
#========================    
    #=============================================
    #This set of plot commands will produce 1 graph with 4 sets of data.
    #Red=Infected Bugs, Pink= Not Infected Bugs, Blue= Infected Cuy
    #Light Blue= Not infected cuy
    #This should be used with only one simulation (nsims=1)
    #============================================= 
    
    #plot(cuyTime, pch="-", col="blue", ylim=c(0, Nbug), main="Red-Inf Bugs, Pink-Not Inf Bugs, Blue-Inf Cuy, Light Blue- Not Inf Cuy")
    #points((Ncuy-cuyTime), pch="-", col="skyblue")
    #points(bugTime, cex=.5, col="red")
    #points((Nbug-bugTime), cex=.5, col="pink")
    
    #=============================================
    #This set of plot commands will produce 4 histograms.  
    #This should be used with many simulations
    #=============================================
                       #Makes the graph split into 4 parts
    setwd(outdir)							#set to the directory to save files
    #setwd("C:/Users/Kate/Desktop/Bug-Cuy Code/No Bottle Analysis 02-11-11")
    #jpeg(paste("No Bottleneck, Run Number",VrunNumber[p],".jpeg"), width = 480, 
    #    height = 480, units = "px", pointsize = 12, quality = 75, bg = "white", res = NA, restoreConsole = TRUE)
   
	jpeg(paste("Test_feb13 Run Number",VrunNumber[p],".jpeg", sep=""), width = 480, height = 480, units = "px", pointsize = 12, quality = 75, bg = "white", res = NA, type="quartz")   
	     
    par(mfrow=c(2,1))
    color<-rainbow(nsims)
    plot(outVector[1,], ylim=c(0, Nbug), main="Number of Vectors Infected Each Day", cex=.2, pch="-")
    for(phi in 2:nsims)
      {
      points(outVector[phi,], cex=.2, pch="-", col=color[phi])
      }
    plot(outHost[1,], ylim=c(0, Ncuy), main="Number of Hosts Infected Each Day", cex=.2, pch="-")
      for(phi in 2:nsims)
      {
      points(outHost[phi,], cex=.2, pch="-", col=color[phi])
      }
     dev.off()
     colnames(outHost)<-paste(c("Day"), 1:EndTime, sep="")
     rownames(outHost)<-paste(c("Simulation"), 1:nsims, sep="")
     colnames(outVector)<-paste(c("Day"), 1:EndTime, sep="")
     rownames(outVector)<-paste(c("Simulation"), 1:nsims, sep="")
     write.csv(outVector, file=paste("Vector Output Run Number",VrunNumber[p],".csv"))
     write.csv(outHost, file=paste("Host Output Run Number",VrunNumber[p],".csv"))
    #hist(out[,1])                          #Top left if out[,1] which is bugs
    #hist(out[,2])                          #Top right is out[,2] cuy
    #hist(bIDmatrix, breaks=50)             #Lower left is bug infect day 
    #hist(cIDmatrix, breaks=50)             #Lower right is cuy infect day 
    
    #==============================================
    #This set of commans automatically saves graphs produced
    #==============================================
    
#    jpeg(paste("No Bottleneck, Run Number ", VrunNumber[p], ".jpeg"), width = 480, height = 480, units = "px", pointsize = 12, 
#        quality = 75, bg = "white", res = NA, restoreConsole = TRUE)
#    par(mfrow=c(2,1))   
#    hist(maxcuy, main=paste("No Bottleneck, Run ", VrunNumber[p]), ylim=c(0,100), sub="Histogram of the max number of infected bugs",
#          breaks=50, xlim=c(0, Ncuy))
#    hist(maxbugs, main="", ylim=c(0,100), sub="Histogram of the max number of infected bugs", 
#          breaks=50, xlim=c(0, Nbug))
#    dev.off()
          }   #Closes for-loop that goes through the runs from the csv


